Self-organized critical neural networks 
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A mechanism for self-organization of the degree of connectivity in model neural networks is studied. 
Network connectivity is regulated locally on the basis of an order parameter of the global dynam- 
ics which is estimated from an observable at the single synapse level. This principle is studied in 
a two-dimensional neural network with randomly wired asymmetric weights. In this class of net- 
works, network connectivity is closely related to a phase transition between ordered and disordered 
dynamics. A slow topology change is imposed on the network through a local rewiring rule moti- 
vated by activity-dependent synaptic development: Neighbor neurons whose activity is correlated, 
on average develop a new connection while uncorrelated neighbors tend to disconnect. As a result, 
robust self-organization of the network towards the order disorder transition occurs. Convergence 
is independent of initial conditions, robust against thermal noise, and does not require fine tuning 
of parameters. 
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Information processing in living organisms is often per- 
formed by large networks of interacting cells with an over- 
all stunning degree of complexity. How can such networks 
be efficiently constructed and how can a robust func- 
tioning be ensured? The observed complexity of many 
nervous systems exceeds by far what can be hard coded 
in the genome [1]. Therefore, developmental principles 
play a key role in network construction. Furthermore, 
as learning is a major function of such networks, self- 
organization and adaptation processes continue through- 
out the lifetime of a network. 

But how can robustness of large dynamical networks 
be ensured in the face of continuous developmental and 
adaptive processes? In general, dynamical stability of 
large networks of dynamical elements and robustness 
against perturbations are not obtained for free: Model 
networks with asymmetric connectivity patterns often 
exhibit regimes of chaotic dynamics with large parameter 
ranges where network dynamics is not easily controlled 
[2]. In networks whose central function is information 
transfer, these regimes would instantly render them use- 
less. Consider, for example, model neural networks with 
asymmetric synaptic couplings, where a percolation tran- 
sition between regimes of ordered and disordered dynam- 
ics is known [3]. In the disordered phase, which occurs 
for densely connected networks, already small perturba- 
tions percolate through the networks.* In such networks, 
developmental processes that change connectivity always 
face the risk of driving the network into the highly con- 
nected regime (where chaotic dynamics prevails), as long 



*This is reminiscent of avalanche like propagation of activity 
in the brain which is observed in some diseases of the central 
nervous system [4]. 



as no explicit mechanism is given that controls the global 
degree of connectivity. 

We here study this question of dynamical robustness 
of networks in the presence of developmental processes in 
the context of a simple toy model, an asymmetric neural 
network combined with simple topology-changing rules. 
In particular we ask how a local rewiring mechanism 
could control global dynamical properties of a large net- 
work and actively contribute to avoiding chaotic regimes. 
While an obvious possibility is a direct feedback of the 
global dynamical state to the synapses, e.g., controlling 
synaptic growth rates, we here consider an even simpler 
mechanism that relies on local information only and in 
principle could be at work in natural systems. We ar- 
gue that if an order parameter characterizing a global 
phase transition is accessible at the single synapse level, 
it can provide the basis for a regulation of global network 
connectivity solely on the basis of local mechanisms. 

Recent models of self-organization of network struc- 
tures show that it is possible to locally measure a global 
order parameter connected to the percolation transition 
of the network, namely the average activity of a single 
node over time [5]. Here we will see that, similarly, the 
average correlation between the activities of two neurons 
contains information about the global order parameter as 
well. The network can then use this approximate order 
parameter to guide the developmental rule. An inter- 
esting question is whether self-organization to a critical 
dynamical transition could occur in a model neural net- 
work on the basis of such a correlation. A possible rule 
is that new synaptic connections preferentially grow be- 
tween correlated neurons, as suggested by the early ideas 
of Hebb [6] and the observation of activity-dependent 
neural development [7]. In the remainder of this paper 
let us study this problem in the framework of a specific 
toy model. We will first define a neural network model 
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with a simple mechanism of synaptic development. Then, 
with numerical studies we will discuss the interplay of 
dynamics on the network with dynamics of the network 
topology. Finally, robustness of self-organizing processes 
in this model and possible implications for biological sys- 
tems are discussed. 

Let us consider a two-dimensional neural network with 
random asymmetric weights on the lattice. The neigh- 
borhood of each neuron is chosen as its Moore neighbor- 
hood with eight neighbors. ^ The weights Wij are ran- 
domly drawn from a uniform distribution Wij G [—1, +1] 
and are nonzero between neighbors, only. Note that 
weights Wij are asymmetric, i.e., in general, w^ ^ Wji. 
Within the neighborhood of a node, a fraction of its 
weights Wij may be set to 0. The average number of 
nonzero weights per node is called the average connec- 
tivity K of the network (for definiteness count e.g. the 
incoming weights at each node, only). The network con- 
sists of N neurons with states a% = ±1 which are updated 
in parallel with a stochastic Little dynamics on the ba- 
sis of inputs received from the neighbor neurons at the 
previous time step: 

prob[<7i(t+l) = +1] =g p (Mt)) 

probMt + 1) = -1] = 1 - g p (fi(t)) (1) 

with 

N 

fi{t) = Y J w i j<j j {t) + e i (2) 

and 

9f>(fiV)) = i + e -2 Ph (t) ( 3 ) 

with the inverse temperature f3 and a threshold 0$. The 
threshold is chosen here as 0i = — 0.1 + 7 and includes 
a small random noise term 7 from a Gaussian of width 
e. This noise term is motivated by the slow fluctuations 
observed in biological neural systems [8] . With respect to 
varying either 9 or K, the network exhibits a percolation 
transition between a phase of ordered dynamics, with 
short transients and short limit cycle attractors, and a 
phase of chaotic dynamics where the length of dynamical 
patterns scales exponentially with system size [3,9]. 



The second part of the model is a slow change of the 
topology of the network by local rewiring of synaptic 
weights: If the activity of two neighbor neurons is on 
average highly correlated (or anticorrclated) , they will 
obtain a common link. If their activity on average is less 
correlated, they will lose their common link. To be more 
specific, let us define the average correlation Cy(r) of a 
pair of neurons over a time interval r 

C ij (T) = -±--f^a i (t)* j (t). (4) 

t=t 

The full model dynamics is then defined as follows. 

1. Start with a random network with an average con- 
nectivity (number of nonzero weights per neuron) 
K ini and a random initial state vector c?(0) = 
(cti(0),...,ctjv(0)). 

2. For each neuron i, choose a random threshold 0i 
from a Gaussian distribution of width e and mean 
fj,. 

3. Starting from the initial state, calculate the new 
system state applying eq. (1) using parallel update. 
Iterate this for r time steps. 

4. Randomly choose one neuron i and one of its 
neighbors j and determine the average correlation 
Cij(r/2) over the last r/2 time steps. (Alter- 
natively, the correlation can be obtained from a 
synaptic variable providing a moving average at any 
given time). 

5. If |Cij(r)| is larger than a given threshold a, i re- 
ceives a new link w^ from site j with a weight cho- 
sen randomly from the interval w^ E [—1,1]} If 
|Cij- (t) I < a, the link w^ is set to (if nonzero). 

6. Go to step 2 and iterate, using the current state of 
the network as new initial state. 

The dynamics of this network is continuous in time, with 
neuron update on a fast time scale and topology update 
of the weights on a well-separated slow "synaptic plas- 
ticity" time scale. Note that the topology-changing rule 
does not involve any global knowledge, e.g., about at- 
tractors. A typical scenario of this dynamical evolution 
is shown in Fig. 1 



^The choice of the type of neighborhood is not critical, how- 
ever, here the Moore neighborhood is more convenient than 

the von Neumann type since, in the latter case, the criti- *Also binary weights could be used as in Ref. [5]. 
cal link density (fraction of nonzero weights) at the percola- 
tion threshold accidentally coincides with the attractor of the 
trivial developmental rule of producing a link with p = 0.5. 
In general, also random sparse neighborhoods would work as 
demonstrated in Ref. [5]. 
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large degree of noise. Figure 3 shows the evolved average 
connectivity as a function of the inverse temperature (3. 
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FIG. 1. Evolution of the average number of nonzero 
weights per neuron over evolutionary time, for a system of 
size N — 64 (8 x 8) and two different initial connectivities 
{Kini = 1.0 and Ki ni = 4.0). Independent of the initial con- 
ditions the networks evolve to a specific average connectivity. 
Parameters are (3 = 25, e = 0.1, a correlation cutoff a = 0.8, 
and an averaging time window of r = 200. 

where the average number of nonzero weights per neu- 
ron K ev is shown as a time series and as cumulative mean. 
One observes that the continuous network dynamics, in- 
cluding the slow local change of the topology, results in 
a convergence of the average connectivity of the network 
to a characteristic value which is independent of initial 
conditions. 
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FIG. 2. Finite size scaling of the evolved average connec- 
tivity. Averages are taken over 4 x 10 5 time steps. 

Finite size scaling of the resulting average connectivity 
indicates the convergence towards a characteristic value 
for large network size iV and exhibits the scaling rela- 
tionship 



K ev (N) = aN- 5 + b 



(5) 



with a = 1.2 ± 0.4, 5 = 0.86 ± 0.07, and b = 2.24 ± 0.03. 
Thus, in the large system size limit N — > oo the net- 
works evolve towards = 2.24 ± 0.03 (see Fig. 2). 
The self-organization towards a specific average connec- 
tivity is largely insensitive to thermal noise of the network 
dynamics, up to ~ 10% of thermal switching errors (or 
(3 > 10) of the neurons. This indicates that the struc- 
ture of a given dynamical attractor is robust against a 




FIG. 3. Evolved average connectivity K ev as a function of 
the inverse temperature /3. Each point is averaged over 10 5 
time steps in a network of size N = 64 and a = 0.5. 

While the stability of dynamical attractors on an in- 
termediate time scale is an important requirement for 
the local sampling of neural correlation, on the long time 
scale of global topological changes, switching between at- 
tractors is necessary to ensure ergodicity at the attractor 
sampling level. The second source of noise, the slow ran- 
dom change in neural thresholds as defined in step (2) of 
the algorithm, is closely related to such transitions be- 
tween attractors. While, in general, the model converges 
also when choosing some arbitrary fixed threshold 6 and 
omitting step (2) from the algorithm, a small threshold 
noise facilitates transitions between limit cycle attractors 
[10] and thus improves sampling over all attractors of a 
network, resulting in an overall increased speed and ro- 
bustness of the convergence. An asynchronous change of 
the threshold 9i, updating one random 9i after complet- 
ing one sweep (time step) of the network, leads to similar 
results as the parallel rule defined above. 

The basic mechanism of the observed self-organization 
in this system is the weak coupling of topological change 
to an order parameter of the global dynamical state of 
the network, and thus is different from the mechanism of 
extremal dynamics, underlying many prominent models 
of self-organized criticality [11]. To illustrate this, let us 
for a moment consider the absolute average correlation 
|CVf( r )l of two neurons which is the parameter used as a 
criterion for the rewiring process. For random networks, 
this quantity is shown in Fig. 4 for different connectivities 
K. 
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FIG. 4. The average correlation ^CV, (t) | between random 
neurons of random networks at different connectivities K. 
Samples are taken over 1000 random networks with 100 ran- 
dom initial conditions each, for network size N = 64. 

Note that the correlation is large for networks with 
small connectivity, and small for networks that are 
densely connected. The rewiring rule balances between 
these two regimes: For high correlation, it is more likely 
that a link is created, at low correlation, links are van- 
ishing. The balance is reached most likely in the region 
of the curve where the slope reaches its maximum, as 
here the observed correlation reacts most sensitively to 
connectivity changes. As the steep portion of the cor- 
relation curve occurs in a region of small connectivities 
where also the critical connectivity K c w 2 of the network 
is located, this makes the correlation measure sensitive 
to the global dynamical state of the network and poten- 
tially useful as an approximation of the order parameter. 
Synaptic development dependent on averaged correlation 
between neurons can thus obtain approximate informa- 
tion about the global dynamical state of the network as is 
realized in the above toy model with a simple implemen- 
tation on the basis of a threshold a. The exact choice of 
the threshold a is not critical, which can be seen from 
the histogram of the absolute correlation |Cy(Y) shown 
in Fig. 5 for a typical run of the model. 



1.2 
1 

0.8 
0.6 
0.4 
0.2 h 




ITTTitTTTTTfl IrmrTlimr^^TrrTTTTt^miTfr 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

FIG. 5. Histogram of | Cij(r) | for a network evolving in 
time, with N = 64 and f3 = 10, taken over a run of 4 x 10 5 
time steps. 



Correlations appear to cluster near high and near low 
values such that the cutoff can be placed anywhere inbe- 
tween the two regimes. Even a threshold value close to 1, 
as compared with the correlation cutoff a = 0.8 used in 
the simulations here, only leads to a minor shift in K ev 
and does not change the overall behavior. 

Up to now we focused on changes of the network struc- 
ture as a result of the dynamics on the network. A further 
aspect is how the structural changes affect the dynam- 
ics on the network itself. Do also dynamical observables 
of the networks self-organize as a result of the observed 
convergence of the network structure? An interesting 
quantity in this respect is the average length of periodic 
attractors as shown in Fig. 6. 
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FIG. 6. Evolution of the cumulative average of attractor 
length for the same system as shown in Fig. 1. The mean 
attractor length converges to a value independent of the two 
initial conditions of the network shown here. The attractor 
length is measured at zero temperature in order to have an 
exactly defined measure. 

Indeed, this dynamical observable of the network dy- 
namics converges to a specific value independent of the 
initial network, similarly to the convergence of the struc- 
tural parameter K considered earlier. From the K depen- 
dency of the neural pair correlation we have seen above 
that the rewiring criterion tends to favor connectivities 
near the critical connectivity of the network. Does also 
the evolved average attractor length relate to critical 
properties of the percolation transition? An approximate 
measure of this aspect is the finite size scaling of the 
evolved average period as shown in Fig. 7. 
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FIG. 7. Finite size scaling of the evolved average attractor 
period (b). Also shown for comparison is the correspond- 
ing scaling of the attractor lengths of an overcritical ran- 
dom network (a) with K = 3.8 and an undercritical one 
(c) with K = 1.5. Symbols denote measured values and 
lines correspond to the fits f a (x) = 15.1a; ' 57 (e 0M5x - l), 
f b (x) = 0.6x 15 , and f c (x) = 0.28a; ' 75 ln(0.097a-). 

For static networks we find that the attractor lengths 
typically scale exponentially with TV in the overcritical 
regime, but less than linearly in the ordered regime. For 
the evolved connectivity K ev in our model, we observe 
scaling close to criticality. Large evolved networks ex- 
hibit relatively short attractors, which otherwise for ran- 
dom networks in the overcritical regime could only be 
achieved by fine tuning. The self-organizing model stud- 
ied here evolves nonchaotic networks without the need 
for parameter tuning. 

To summarize, neural network development has been 
studied in an asymmetric model neural network. The de- 
velopmental rule is based on local rewiring motivated by 
Hebbian, activity-dependent synaptic development. In 
a continuously running network, robust self-organization 
of the network towards the percolation transition be- 
tween ordered and disordered dynamics is observed, in- 
dependent of initial conditions and robust against ther- 
mal noise. The basic model is robust against changes in 
the details of the algorithm. We conclude that a weak 
coupling of the rewiring process to an approximate mea- 
surement of an order parameter of the global dynamics 
is sufficient for a robust self-organization towards criti- 
cality. In particular, the order parameter has been es- 
timated solely from information available on the single 
synapse level via time averaging of correlated neural ac- 
tivities. 

While here we considered self-organization in model 
neural networks, the observed mechanism may occur in 
other more complex systems. For example, global dy- 
namical order from self-organization at the synapse level 
could, in principle, be at work in biological nervous sys- 
tems as well. Prerequisites are an averaging procedure 
of correlated activities on slow time scales (similar to 
synaptic processes underlying learning through long term 
potentiation), and a coupling to synaptic development. 
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